﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Triangulator.Geometry;
using Extreme.Mathematics.LinearAlgebra;
using Extreme.Mathematics;


namespace erse253_homework
{
    public class PointEstimation
    {

        private double EstimateByOrdinaryKrigingMethod(List<PointxyV> vectorData, Matrix invertedMatrix, int i, int j)
        {

            Vector DMatrix = Vector.Create(vectorData.Count + 1);
            for (int k = 0; k < vectorData.Count; k++)
                DMatrix[k] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[k].i - i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[k].j - j) / 10, 2) / 27.0);

            DMatrix[vectorData.Count] = 1;

            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);


            double estimatedValue = 0;

            for (int s = 0; s < omega.Length - 1; s++)
                estimatedValue += vectorData[s].V * omega[s];

            return estimatedValue;

        }

        private double EstimateBySimpleKrigingMethod(List<PointxyV> vectorData, Matrix invertedMatrix, int i, int j)
        {

            Vector DMatrix = Vector.Create(vectorData.Count);
            for (int k = 0; k < vectorData.Count; k++)
                DMatrix[k] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[k].i - i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[k].j - j) / 10, 2) / 27.0);


            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);


            double estimatedValue = 0;

            for (int s = 0; s < omega.Length; s++)
                estimatedValue += vectorData[s].V * omega[s];

            return estimatedValue;

        }

        private double EstimateValuePolygonMethod(List<PointxyV> points, int i, int j)
        {
            for (int l = 0; l < points.Count; l++)
                points[l].UpdateDistance(i, j);

            points = OrderPoints(points);
            return points[0].V;
            
        }

        private double EstimateValueByInverseSquaredDistance(List<PointxyV> points, int i, int j)
        {
         
                double normalize = 0;
                double estimatedValue = 0;
                double distance = 0;
                for (int l = 0; l < points.Count; l++)
                {
                    points[l].UpdateDistance(i, j);
                    distance = 1.0 / (Math.Pow(points[l].distance, 2) / 100);
                    estimatedValue += points[l].V*distance;
                    normalize += distance;

                }

                return estimatedValue / normalize;

        }

        private double EstimateValueTriangulationMethod(double[,] matrix, List<Point> Vertex, List<Triangle> triangles, int i, int j)
        {

                foreach (Triangulator.Geometry.Triangle triangle in triangles)
                {
                    List<PointxyV> points = SetDelaunyTriangle(matrix, Vertex, i, j, triangle);

                    if (PointInTriangle(points, new PointxyV(i, j, matrix[i, j], 0, 0)))
                    {
                        Matrix dataMatrix = Matrix.Create(LoadMatrix(points));
                            
                        Vector dataVector = Vector.Create(LoadArray(points));

                        Vector answer = dataMatrix.Solve(dataVector);
                        return answer[0] * ((double)i) / 10 + answer[1] * ((double)j) / 10 + answer[2];
                    }
                }
                return double.NaN;
        }

        private double EstimateErrorVarianceOrdinaryKrigingMethod(List<PointxyV> vectorData, Matrix invertedMatrix, int i, int j)
        {

            Vector DMatrix = Vector.Create(vectorData.Count + 1);
            for (int k = 0; k < vectorData.Count; k++)
                DMatrix[k] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[k].i - i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[k].j - j) / 10, 2) / 27.0);

            DMatrix[vectorData.Count] = 1;

            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);


            double estimatedValue = 0;

            for (int s = 0; s < omega.Length - 1; s++)
                estimatedValue += omega[s] * DMatrix[s];


            //must return covariance of the element minus estimated value plus mi.
            return 32.0 - estimatedValue + omega[omega.Length - 1];

        }

        private double EstimateErrorVarianceSimpleKrigingMethod(List<PointxyV> vectorData, Matrix invertedMatrix, int i, int j)
        {

            Vector DMatrix = Vector.Create(vectorData.Count);
            for (int k = 0; k < vectorData.Count; k++)
                DMatrix[k] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[k].i - i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[k].j - j) / 10, 2) / 27.0);

            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);


            double estimatedValue = 0;

            for (int s = 0; s < omega.Length; s++)
                estimatedValue += omega[s] * DMatrix[s];


            //must return covariance of the element minus estimated value.
            return 32.0 - estimatedValue;

        }


        #region Cross Validation
        public void CrossValidateKriging(double[,] matrix, double[,] estimatedValues)
        {
            Matrix coefficients;

            Parallel.For(0, estimatedValues.GetLength(0), delegate(int i)
            {
                List<PointxyV> vectorData;
                PointxyV temp = new PointxyV();

                GetDataPoints(matrix, out vectorData);


                for (int j = 0; j < estimatedValues.GetLength(1); j++)
                {
                    if (!double.IsNaN(matrix[i, j]))
                    {


                        for (int k = 0; k < vectorData.Count; k++)
                            if (vectorData[k].i == i && vectorData[k].j == j)
                            {
                                temp = vectorData[k];
                                vectorData.RemoveAt(k);
                                break;
                            }

                        coefficients = Matrix.Create(CreateCovarianceMatrix(vectorData));

                        Matrix invertedMatrix = Matrix.Inverse(coefficients);

                        estimatedValues[i, j] = EstimateByOrdinaryKrigingMethod(vectorData, invertedMatrix, i, j);
                        vectorData.Add(temp);
                    }
                }
            });

        }

        public void CrossValidatePolygon(double[,] matrix, double[,] estimatedValues)
        {
            System.Threading.Tasks.Parallel.For(0, matrix.GetLength(0), delegate(int i)
            {
                List<PointxyV> points = GetPointList(estimatedValues);
                PointxyV temp = new PointxyV();
                for (int j = 0; j < matrix.GetLength(1); j++)
                    if (!double.IsNaN(matrix[i, j]))
                    {
                        //revome the point i j from list then add again after the function executes.
                        for (int k = 0; k < points.Count; k++)
                            if (points[k].i == i && points[k].j == j)
                            {
                                temp = points[k];
                                points.RemoveAt(k);
                                break;
                            }
                        estimatedValues[i, j] = EstimateValuePolygonMethod(points, i, j);
                        points.Add(temp);
                    }
            });
        }

        public void CrossValidateTriangulation(double[,] matrix, double[,] estimatedValues)
        {
            System.Threading.Tasks.Parallel.For(0, matrix.GetLength(0), delegate(int i)
            {

                for (int j = 0; j < matrix.GetLength(1); j++)
                    if (!double.IsNaN(matrix[i, j]))
                    {
                        double[,] estimateElement = GetMatrixInUse(matrix);
                        estimateElement[i, j] = double.NaN;


                        DelaunyTriangulation delauny = new DelaunyTriangulation();
                        List<Point> Vertex = new List<Point>();

                        for (int k = 0; k < estimateElement.GetLength(0); k++)
                            for (int l = 0; l < estimateElement.GetLength(1); l++)
                                if (!double.IsNaN(estimateElement[k, l]))
                                    Vertex.Add(new Point(k, l));

                        List<Triangle> triangles = delauny.Triangulate(Vertex);

                        estimatedValues[i, j] = EstimateValueTriangulationMethod(estimateElement, Vertex, triangles, i, j);
                    }
            });
        }

        public void CrossValidateInverseSquaredDistance(double[,] matrix, double[,] estimatedValues)
        {
            Parallel.For(0, estimatedValues.GetLength(0), delegate(int i)
            {
                List<PointxyV> points = GetPointList(matrix);
                PointxyV temp = new PointxyV();
                for (int j = 0; j < matrix.GetLength(1); j++)
                    if (!double.IsNaN(matrix[i, j]))
                    {

                        for (int k = 0; k < points.Count; k++)
                            if (points[k].i == i && points[k].j == j)
                            {
                                temp = points[k];
                                points.RemoveAt(k);
                                break;
                            }

                        estimatedValues[i, j] = EstimateValueByInverseSquaredDistance(points, i, j);
                        points.Add(temp);

                    }
            });
        }
        #endregion

        #region AuxiliaryMethods

        private void SetLagrangeCondidition(double[,] coefficients)
        {

            for (int j = 0; j < coefficients.GetLength(1); j++)
                coefficients[coefficients.GetLength(0) - 1, j] = 1;


            for (int i = 0; i < coefficients.GetLength(0); i++)
                coefficients[i, coefficients.GetLength(1) - 1] = 1;


            coefficients[coefficients.GetLength(0) - 1, coefficients.GetLength(1) - 1] = 0;
        }

        private double[,] LoadMatrix(List<PointxyV> triangle)
        {
            double[,] matrix = new double[3, 3];

            for (int i = 0; i < 3; i++)
            {
                matrix[i, 0] = ((double)triangle[i].i) / 10;
                matrix[i, 1] = ((double)triangle[i].j) / 10;
                matrix[i, 2] = 1;
            }

            return matrix;
        }

        private double[] LoadArray(List<PointxyV> triangle)
        {
            double[] vector = new double[3];

            for (int i = 0; i < 3; i++)
            {
                vector[i] = triangle[i].V;
            }

            return vector;
        }

        private float sign(PointxyV p1, PointxyV p2, PointxyV p3)
        {
            return (p1.i - p3.i) * (p2.j - p3.j) - (p2.i - p3.i) * (p1.j - p3.j);
        }

        private bool PointInTriangle(List<PointxyV> p, PointxyV pt)
        {
            bool b1, b2, b3;
            //todo here!
            b1 = sign(pt, p[0], p[1]) <= 0.0f;
            b2 = sign(pt, p[1], p[2]) <= 0.0f;
            b3 = sign(pt, p[2], p[0]) <= 0.0f;

            return ((b1 == b2) && (b2 == b3));
        }

        private double[,] GetMatrixInUse(double[,] m)
        {
            double[,] matrix = new double[m.GetLength(0), m.GetLength(1)];
            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = 0; j < matrix.GetLength(1); j++)
                    matrix[i, j] = m[i, j];

            return matrix;

        }

        private List<PointxyV> CreatePointList(PointxyV p, PointxyV p2, PointxyV p3)
        {
            List<PointxyV> points2 = new List<PointxyV>();
            points2.Clear();
            points2.Add(p);
            points2.Add(p2);
            points2.Add(p3);

            return points2;
        }

        private List<PointxyV> GetPointList(int tolerance, int i, int j, double[,] matrix)
        {

            List<PointxyV> points = new List<PointxyV>();

            for (int k = -tolerance / 2; k < tolerance / 2; k++)
                for (int l = -tolerance / 2; l < tolerance / 2; l++)
                    if ((i + k >= 0) && (j + l >= 0) && (i + k < matrix.GetLength(0)) && (j + l < matrix.GetLength(1)) && !double.IsNaN(matrix[i + k, j + l]))
                        points.Add(new PointxyV(i + k, j + l, matrix[i + k, j + l], i, j));


            return points;
        }

        private List<PointxyV> GetPointList(double[,] matrix)
        {

            List<PointxyV> points = new List<PointxyV>();

            for (int k = 0; k < matrix.GetLength(0); k++)
                for (int l = 0; l < matrix.GetLength(1); l++)
                    if (!double.IsNaN(matrix[k, l]))
                        points.Add(new PointxyV(k, l, matrix[k, l], 0, 0));


            return points;
        }

        private double[,] CreateCovarianceMatrix(List<PointxyV> vectorData)
        {
            double[,] coefficients = new double[vectorData.Count + 1, vectorData.Count + 1];

            System.Threading.Tasks.Parallel.For(0, vectorData.Count, delegate(int i)
            {
                for (int j = 0; j < vectorData.Count; j++)
                    coefficients[i, j] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[i].i - vectorData[j].i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[i].j - vectorData[j].j) / 10, 2) / 27.0);

            });


            //improve condition number by adding 1 to the diagonal elements.
            for (int i = 0; i < coefficients.GetLength(0) - 1; i++)
                for (int j = 0; j < coefficients.GetLength(1) - 1; j++)
                {
                    if (i == j)
                        coefficients[i, j] += 1;
                }

            SetLagrangeCondidition(coefficients);
            return coefficients;
        }

        private double[,] CreateCovarianceMatrixSimpleKriging(List<PointxyV> vectorData)
        {
            double[,] coefficients = new double[vectorData.Count, vectorData.Count];

            System.Threading.Tasks.Parallel.For(0, vectorData.Count, delegate(int i)
            {
                for (int j = 0; j < vectorData.Count; j++)
                    coefficients[i, j] = 32.0 * Math.Exp(-Math.Pow((double)(vectorData[i].i - vectorData[j].i) / 10, 2) / 35.0 - Math.Pow((double)(vectorData[i].j - vectorData[j].j) / 10, 2) / 27.0);

            });


            //improve condition number by adding 1 to the diagonal elements.
            for (int i = 0; i < coefficients.GetLength(0); i++)
                for (int j = 0; j < coefficients.GetLength(1); j++)
                {
                    if (i == j)
                        coefficients[i, j] += 1;
                }

            return coefficients;
        }

        public void GetDataPoints(double[,] matrix, out List<PointxyV> vectorData)
        {
            vectorData = new List<PointxyV>();
            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = 0; j < matrix.GetLength(1); j++)
                    if (!double.IsNaN(matrix[i, j]))
                        vectorData.Add(new PointxyV(i, j, matrix[i, j], 0, 0));


            //vectorData = OrderPoints(vectorData);
        }

        private static List<PointxyV> SetDelaunyTriangle(double[,] matrix, List<Point> Vertex, int i, int j, Triangle triangle)
        {
            List<PointxyV> points = new List<PointxyV>();
            points.Add(new PointxyV((int)Vertex[triangle.p1].X, (int)Vertex[triangle.p1].Y, matrix[(int)Vertex[triangle.p1].X, (int)Vertex[triangle.p1].Y], i, j));
            points.Add(new PointxyV((int)Vertex[triangle.p2].X, (int)Vertex[triangle.p2].Y, matrix[(int)Vertex[triangle.p2].X, (int)Vertex[triangle.p2].Y], i, j));
            points.Add(new PointxyV((int)Vertex[triangle.p3].X, (int)Vertex[triangle.p3].Y, matrix[(int)Vertex[triangle.p3].X, (int)Vertex[triangle.p3].Y], i, j));

            return points;
        }

        private List<PointxyV> OrderPoints(List<PointxyV> points)
        {
            return points.OrderBy(point => point.distance).ToList();
        }

        public List<PointxyV> GetPointsIntoBlock(double[,] matrix, int i, int j, int sizeX, int sizeY)
        {
            List<PointxyV> pointsInsideArea = new List<PointxyV>();

            for (int k = 0; k < sizeX; k++)
                for (int l = 0; l < sizeY; l++)
                    if (!double.IsNaN(matrix[k + i, l + j]))
                        pointsInsideArea.Add(new PointxyV(k + i, l + j, matrix[k + i, l + j], i, j));

            return pointsInsideArea;
        }

        private List<PointxyV> GetPointsIntoArea(double[,] vectorData, int size, int i, int j)
        {

            List<PointxyV> pointsInsideArea = new List<PointxyV>();

            for (int k = i - size / 2; k < i + size / 2; k++)
                for (int l = j - size / 2; l < j + size / 2; l++)
                    if (k >= 0 && k < vectorData.GetLength(0) && l >= 0 && l < vectorData.GetLength(1) && !double.IsNaN(vectorData[k, l]))
                        pointsInsideArea.Add(new PointxyV(k, l, vectorData[k, l], i, j));

            return pointsInsideArea;
        }
        #endregion

        public double[,] OrdinaryKrigingMethod(double[,] SparceMatrix)
        {
            List<PointxyV> vectorData;
            double[,] estimatedMatrix = SparceMatrix;
            GetDataPoints(estimatedMatrix, out vectorData);


            Matrix coefficients = Matrix.Create(CreateCovarianceMatrix(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {
                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    if (double.IsNaN(estimatedMatrix[i, j]))
                        estimatedMatrix[i, j] = EstimateByOrdinaryKrigingMethod(vectorData, invertedMatrix, i, j);
            });
            return estimatedMatrix;
        }

        public double[,] TrianglulationMethod(double[,] SparceMatrix, double[,] copy)
        {
            double[,] matrix = SparceMatrix;
            double[,] estimatedValues = copy;

            DelaunyTriangulation delauny = new DelaunyTriangulation();
            List<Triangulator.Geometry.Point> Vertex = new List<Triangulator.Geometry.Point>();

            for (int i = 0; i < matrix.GetLength(0); i++)
                for (int j = 0; j < matrix.GetLength(1); j++)
                    if (!double.IsNaN(matrix[i, j]))
                        Vertex.Add(new Triangulator.Geometry.Point(i, j));

            List<Triangulator.Geometry.Triangle> triangles = delauny.Triangulate(Vertex);
            PointEstimation triangulation = new PointEstimation();
            System.Threading.Tasks.Parallel.For(0, estimatedValues.GetLength(0), delegate(int i)
            {
                for (int j = 0; j < estimatedValues.GetLength(1); j++)
                    if (double.IsNaN(estimatedValues[i, j]))
                        estimatedValues[i, j] = triangulation.EstimateValueTriangulationMethod(matrix, Vertex, triangles, i, j);
            });
            return estimatedValues;
        }

        public double[,] InverseSquaredDistance(double[,] SparceMatrix, double[,] copy)
        {
            double[,] estimatedData = SparceMatrix;
            double[,] dataPoints = copy;
            PointEstimation inverseSquaredDistance = new PointEstimation();

            Parallel.For(0, dataPoints.GetLength(0), delegate(int i)
            {
                List<PointxyV> points = GetPointList(dataPoints);
                for (int j = 0; j < dataPoints.GetLength(1); j++)
                    if (double.IsNaN(dataPoints[i, j]))
                        estimatedData[i, j] = inverseSquaredDistance.EstimateValueByInverseSquaredDistance(points, i, j);
            });
            return estimatedData;
        }

        public double[,] PolygonMethod(double[,] SparceMatrix, double[,] copy)
        {
            //get the data matrix and a copy of it to calculate intermediate resutls.
            double[,] matrix = SparceMatrix;
            double[,] m = copy;
            PointEstimation polygon = new PointEstimation();

            System.Threading.Tasks.Parallel.For(0, matrix.GetLength(0), delegate(int i)
            {
                List<PointxyV> points = GetPointList(m);
                for (int j = 0; j < matrix.GetLength(1); j++)
                    matrix[i, j] = polygon.EstimateValuePolygonMethod(points, i, j);
            });
            return matrix;
        }

        public double[,] SimpleKrigingMethod(double[,] SparceMatrix)
        {
            PointEstimation simpleKriging = new PointEstimation();
            List<PointxyV> vectorData;
            double[,] estimatedMatrix = SparceMatrix;
            simpleKriging.GetDataPoints(estimatedMatrix, out vectorData);


            Matrix coefficients = Matrix.Create(simpleKriging.CreateCovarianceMatrixSimpleKriging(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {
                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    if (double.IsNaN(estimatedMatrix[i, j]))
                        estimatedMatrix[i, j] = simpleKriging.EstimateBySimpleKrigingMethod(vectorData, invertedMatrix, i, j);
            });
            return estimatedMatrix;
        }

        public double EstimateVarianceKriging( List<PointxyV> vectorData)
        {
            Matrix coefficients = Matrix.Create(CreateCovarianceMatrix(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Vector DMatrix = Vector.Create(vectorData.Count + 1);
            for (int i = 0; i < DMatrix.Length; i++)
                DMatrix[i] = 0;

            DMatrix[vectorData.Count] = 1;


            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);

            double estimatedMean = 0;
            for (int s = 0; s < omega.Length - 1; s++)
                estimatedMean += vectorData[s].V * omega[s];

            double estimatedVar = 0;
            for (int s = 0; s < omega.Length - 1; s++)
                estimatedVar += Math.Pow((vectorData[s].V - estimatedMean), 2) * omega[s];
            return estimatedVar;
        }

        public  double EstimateMeanKriging( List<PointxyV> vectorData)
        {

            Matrix coefficients = Matrix.Create(CreateCovarianceMatrix(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Vector DMatrix = Vector.Create(vectorData.Count + 1);
            for (int i = 0; i < DMatrix.Length; i++)
                DMatrix[i] = 0;

            DMatrix[vectorData.Count] = 1;


            Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);

            double estimatedValue = 0;
            for (int s = 0; s < omega.Length - 1; s++)
                estimatedValue += vectorData[s].V * omega[s];

            return estimatedValue;
        }

        public double EstimateMeanUsingBlockOrdinaryKriging(double[,] estimatedMatrix, double[,] dataMatrix, int numberDivisions)
        {
            PointEstimation kriging = new PointEstimation();

            int sizeX = estimatedMatrix.GetLength(0) / numberDivisions;
            int sizeY = estimatedMatrix.GetLength(1) / numberDivisions;

            double blockKrigingMean = 0;
            int number = 0;
            //Parallel.For(0, numberDivisions, delegate(int i){
            for (int i = 0; i < numberDivisions; i++)
                for (int j = 0; j < numberDivisions; j++)
                {
                    List<PointxyV> pointsInsideArea = GetPointsIntoBlock(dataMatrix, i * sizeX, j * sizeY, sizeX, sizeY);

                    if (pointsInsideArea.Count > 0)
                    {
                        Matrix coefficients = Matrix.Create(kriging.CreateCovarianceMatrix(pointsInsideArea));
                        Matrix invertedMatrix = Matrix.Inverse(coefficients);
                        Vector DMatrix = Vector.Create(pointsInsideArea.Count + 1);

                        for (int t = 0; t < DMatrix.Length; t++)
                            DMatrix[t] = 0;

                        DMatrix[pointsInsideArea.Count] = 1;


                        Vector omega = Matrix.Multiply(invertedMatrix, DMatrix);

                        double estimatedValue = 0;
                        for (int s = 0; s < omega.Length - 1; s++)
                            estimatedValue += pointsInsideArea[s].V * omega[s];
                        blockKrigingMean += estimatedValue;
                        number++;

                    }

                }
            blockKrigingMean /= number;
            return blockKrigingMean;
        }

        public double[,] OrdinaryKrigingUsingSearchCriteria(int size, double[,] estimatedMatrix, double[,] dataMatrix)
        {
            PointEstimation kriging = new PointEstimation();

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {

                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    if (double.IsNaN(estimatedMatrix[i, j]))
                    {
                        List<PointxyV> pointsInsideArea = GetPointsIntoArea(dataMatrix, size, i, j);
                        if (pointsInsideArea.Count >= 1)
                        {
                            Matrix coefficients = Matrix.Create(kriging.CreateCovarianceMatrix(pointsInsideArea));
                            Matrix invertedMatrix = Matrix.Inverse(coefficients);
                            estimatedMatrix[i, j] = kriging.EstimateByOrdinaryKrigingMethod(pointsInsideArea, invertedMatrix, i, j);
                        }
                        pointsInsideArea.Clear();
                    }
            });
            return estimatedMatrix;
        }

        public double[,] SimpleKrigingUsingSearchCriteria(int size, double[,] estimatedMatrix, double[,] dataMatrix)
        {

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {

                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    if (double.IsNaN(estimatedMatrix[i, j]))
                    {
                        List<PointxyV> pointsInsideArea = GetPointsIntoArea(dataMatrix, size, i, j);
                        if (pointsInsideArea.Count >= 1)
                        {
                            Matrix coefficients = Matrix.Create(CreateCovarianceMatrixSimpleKriging(pointsInsideArea));
                            Matrix invertedMatrix = Matrix.Inverse(coefficients);
                            estimatedMatrix[i, j] = EstimateBySimpleKrigingMethod(pointsInsideArea, invertedMatrix, i, j);
                        }
                        pointsInsideArea.Clear();
                    }
            });
            return estimatedMatrix;
        }

        public double[,] OrdinaryKrigingGridEstimateErrorVariance(double[,] SparceMatrix)
        {
            List<PointxyV> vectorData;
            double[,] estimatedMatrix = SparceMatrix;
            GetDataPoints(estimatedMatrix, out vectorData);


            Matrix coefficients = Matrix.Create(CreateCovarianceMatrix(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {
                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    estimatedMatrix[i, j] = EstimateErrorVarianceOrdinaryKrigingMethod(vectorData, invertedMatrix, i, j);

            });
            return estimatedMatrix;
        }

        public double[,] SimpleKrigingGridEstimateErrorVariance(double[,] SparceMatrix)
        {
            List<PointxyV> vectorData;
            double[,] estimatedMatrix = SparceMatrix;
            GetDataPoints(estimatedMatrix, out vectorData);


            Matrix coefficients = Matrix.Create(CreateCovarianceMatrixSimpleKriging(vectorData));
            Matrix invertedMatrix = Matrix.Inverse(coefficients);

            Parallel.For(0, estimatedMatrix.GetLength(0), delegate(int i)
            {
                for (int j = 0; j < estimatedMatrix.GetLength(1); j++)
                    estimatedMatrix[i, j] = EstimateErrorVarianceSimpleKrigingMethod(vectorData, invertedMatrix, i, j);

            });
            return estimatedMatrix;
        }

        public void EstimateBlockOrdinaryKriging(double[,] estimatedMatrix, double[,] dataMatrix, PointEstimation kriging, int numberDivisions)
        {
            int sizeX = estimatedMatrix.GetLength(0) / numberDivisions;
            int sizeY = estimatedMatrix.GetLength(1) / numberDivisions;

            //Parallel.For(0, numberDivisions, delegate(int i){
            for (int i = 0; i < numberDivisions; i++)
                for (int j = 0; j < numberDivisions; j++)
                {
                    List<PointxyV> pointsInsideArea = GetPointsIntoBlock(dataMatrix, i * sizeX, j * sizeY, sizeX, sizeY);
                    for (int k = 0; k < sizeX; k++)
                        for (int l = 0; l < sizeY; l++)
                            if (double.IsNaN(estimatedMatrix[i * sizeX + k, j * sizeY + l]) && pointsInsideArea.Count >= 1)
                            {
                                Matrix coefficients = Matrix.Create(kriging.CreateCovarianceMatrix(pointsInsideArea));
                                Matrix invertedMatrix = Matrix.Inverse(coefficients);
                                estimatedMatrix[i * sizeX + k, j * sizeY + l] = kriging.EstimateByOrdinaryKrigingMethod(pointsInsideArea, invertedMatrix, i * sizeX + k, j * sizeY + l);
                            }

                }
            // });
        }

        public void EstimateBlockSimpleKriging(double[,] estimatedMatrix, double[,] dataMatrix, PointEstimation kriging, int numberDivisions)
        {
            int sizeX = estimatedMatrix.GetLength(0) / numberDivisions;
            int sizeY = estimatedMatrix.GetLength(1) / numberDivisions;

            Parallel.For(0, numberDivisions, delegate(int i)
            {

                for (int j = 0; j < numberDivisions; j++)
                {
                    List<PointxyV> pointsInsideArea = GetPointsIntoBlock(dataMatrix, i, j, sizeX, sizeY);
                    for (int k = 0; k < sizeX; k++)
                        for (int l = 0; l < sizeY; l++)
                            if (double.IsNaN(estimatedMatrix[i * sizeX + k, j * sizeY + l]))
                            {

                                if (pointsInsideArea.Count >= 1)
                                {
                                    Matrix coefficients = Matrix.Create(kriging.CreateCovarianceMatrixSimpleKriging(pointsInsideArea));
                                    Matrix invertedMatrix = Matrix.Inverse(coefficients);
                                    estimatedMatrix[i, j] = kriging.EstimateBySimpleKrigingMethod(pointsInsideArea, invertedMatrix, i, j);
                                }
                                pointsInsideArea.Clear();
                            }
                }
            });
        }

    }
}
